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Adsorption of hydrogen atoms on a single graphite sheet (graphene) has been investigated by first- 
principles electronic structure means, employing plane-wave based, periodic density functional 
theory. A reasonably large 5x5 surface unit cell has been employed to study single and multiple 
adsorption of H atoms. Binding and barrier energies for sequential sticking have been computed 
for a number of configurations involving adsorption on top of carbon atoms. We find that binding 
energies per atom range from ~ 0.8 eV to ~ 1.9 eV, with barriers to sticking in the range 0.0 — 0.2 
eV. In addition, depending on the number and location of adsorbed hydrogen atoms, we find 
that magnetic structures may form in which spin density localizes on a \/3x\/3R30° sublattice, 
and that binding (barrier) energies for sequential adsorption increase (decrease) linearly with the 
site-integrated magnetization. These results can be rationalized with the help of the valence-bond 
resonance theory of planar tt conjugated systems, and suggest that preferential sticking due to 
barrierless adsorption is limited to formation of hydrogen pairs. 



I. INTRODUCTION 



Recent years have witnessed an ever growing inter- 
est in carbon-based materials. Carbon, being a small 
atom with a half-filled shell, is able to mix its valence s 
and p orbitals to various degrees, thereby forming the 
building block for extended structures of incredibly 
different electronic, magnetic and mechanical prop- 
erties. Among them, those formed by sp^ C atoms 
have attracted much attention in the last few years. 
They can be collective termed as graphitic compounds 
and comprise graphite, carbon nanotubes, fullerenes, 
Polycyclic Aromatic Hydrocarbons (PAHs), and re- 
cently graphene (the one-atom thick layer of graphite) 
and graphene nanoribbons (GNRs). In particular, the 
revolutionary (and embarrassing simple) fabrication 
of graphene (Novoselov et al. |2004 1 has opened the 



way for a wealth of studies in both fundamental and 
applied science. New, extraordinary properties have 
become available to material design since its isolation. 
Indeed, even though they have been known since th e 
first theoretical analysis by Wallace (Wallace 19471, 



it was only the experimental observation of the ex- 



istence of one-atom tick layer of graphite that trig- 
gered much of the current interest. In particular, one 
of the most interesting aspects of graphene is that 
it presents low energy excitations as massless, chi- 
ral, Dirac fermions mimicking the physics of quantum 



electrodynamics (Novoselov et al. 2005 Zhang et al. 



2005 Castro Neto et al. 2008 1 



In this context, adsorption of hydrogen atoms on 
graphene and GNRs can be used to tailor electronic 
and magnetic properties, as already suggested for 
other 'defects', with the advantage of being much eas- 
ier to realize than e.g. vacancies. In addition, inter- 
action of hydrogen atoms with graphitic compounds 
has been playing an important role in a number of 
fields as diverse as nuclear fusion (Meregalli and Par- 
|rinello| |2001 Mayer et al. 20011, hydrogen storage 



(Schlapbach and Ziittel 



200111 and interstellar chem- 



lams 



19991 



istry ( [Hartquist and Wil 

In material design for hydrogen storage, several car- 
bon based structures has been proposed as candidates 
( Schlapbach and Ziittel 2001 1 , in particular in con- 
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nection with the spillover effect following embedding 
of metallic nanoparticles. Though these materials are 
in practice still far from the weight percent target 
stated by the US department of Energy, they remain a 
cheap and safe alternative, and a deeper understand- 
ing of the mechanisms underlying adsorption may lead 



in future to a more efficient material design. 

In interstellar chemistry hydrogen-graphite and 
hydrogen-PAHs systems have become realistic mod- 
els to investigate molecular hydrogen formation in the 
interstellar medium (ISM). There are still open ques- 
tions in this context since, in spite of continuous de- 
struction by UV radiation and cosmic rays, H2 is the 
most abundant molecule of the ISM. It is now widely 
accepted that H2 can only form on the surface of inter- 



stellar dust grains and particles (Gould and Salpeter 



1963 Hollenbach and Salpeter 1970 19711, which 



with the exception of cold, dense molecular clouds- are 
either carbon-coated silicate grains or carbonaceous 



particles or 


large PAHs ( 


Greenberg 2002 




Williams 


and Herbst 




2002 


Draine 




20031. This fi 


nding has 



stimulated a number of theoretical ( Jeloaica and Sidis 
T999{ ISha and Jackson] [2002] ISha et al.| |2002| |Zecho| 



et al.| [2002] |Sha et al.[ [2005] porisset et al [2QQ4| 



2005[ 



Allouche e t al.| |2005| [Morisset et al 
nazzo and Tantardini , 2005 ; Allouc he et al, 
win et al.l j2006f iMartinazzo and Tantardini, 2006a b 



Marti- 
2006, Ker^ 



Bonfanti et al.[ [200 7 Cuppen and Hornekasr, 2008, 
Medina and Jackson[ [2008 1 and experimental ( [Zechq 



erties. Complementing previous investigations, how- 
ever, we show how the simple tt resonance chemical 
model helps in rationalizing the findings. A parallel 
work on diff erent graphitic substra tes (PAHs) will fol- 
low shortly (Bonfanti et al. 20081. 

The paper is organized as follows. Details of our 
first-principles calculations are given in Section[ll| and 
their results in Section [Illj where we analyze adsorp- 
tion of a single H atom and briefly introduce the chem- 
ical model (Section III. A I, we consider formation of 
pairs (Section |lII.B[ | and formation of three- and four- 
atom clusters (Section III.C I. We summarize and con- 
clude in Section II VI 



II. COMPUTATIONAL METHODS 

Periodic density functional theory as implemented 
in the Vienna Ab initio Simulation Package suite 
(VASP) ([Kresse and Hafner[ [1994[ [1993] [Kresse and" 



Furthmiiller 



et al.[ [2002[ [Giittler etal 



Giittler et al.j [2004b| [Andree et al 



culations. 



2004a[ [Zecho et al.| 120041 



2006 



pornekffirl 



1996a|b I has been used in all the cal 
The projector-augmented wave method 
within the frozen core approximation has been used 
to describe the electron-core interaction ( Blochl[ 1994 



2006^ Islam et al. 2007 Hornekaer et al. 2007| studies 



Kresse and Joubert 19991, with a Perdew-Burke- 



et al., 20 06a|b] [Baouche et al.| [2006[ jCreighan et aTT Ernzerhof (PBE) ([Perdew et al.[ [T^ [ij^ func 



on hydrogen graphitic systems aimed at elucidating 
the possible reaction pathways leading ultimately to 
molecule formation. 

One interesting finding of these studies is the ten- 
dency of hydrogen atoms to cluster at all but very low 



coverage conditions (Andree et al. 2006 Hornekaer 



et al. 2006a|b 20071. New mechanisms for hydrogen 



sticking (Hornek^r et al. 2006b| and new recombi- 
nation pathways (Hornek^r et al. 2006a I have been 



tional within the generalized gradient approximation 
(GGA). Due to the crucial role that spin plays in this 
system all our calculations have been performed in a 
spin unrestricted framework. 

All calculations have used an energy cutoff of 500 
eV and a 6x6x1 F-centered fc-points mesh to span the 
electron density, in a way to include all the special 
points of the cell. The linear tetrahedron method with 
Blochl corrections is used (Blochl et al. 1994 1 together 



proposed, based on the now common agreement that 
the presence of one or more adsorbate atoms strongly 
influences subsequent adsorption. It is clear that such 
an influence can only result as a consequence of a 
substrate-mediated interaction which makes use of the 
unusual electronic properties of graphitic compounds, 
but at present a comprehensive model for multiple 
chemisorption is still missing. 

In this work we present first principles calculations 
of single and multiple adsorption of hydrogen atoms 
on a graphene sheet, used as a model graphitic ma- 
terial, with the aim of understanding the relation- 
ship between the substrate electronic properties and 
the stability of various cluster configurations. This 
work parallels analogous investigations of defects in 
graphene and GNRs (Pereira et al. 2006 Yazyev and 
Helm[ [2007| [Pisani et al.[ [2008[ [Pereira et al.[ [2008| 



with a 0.2 eV smearing. All the atomic positions have 
been fully relaxed until the Hellmann-Feynman forces 
dropped below 10~^ eVA~^, while convergence of the 
electronic structures has been ensured by forcing the 
energy difference in the self consistent cycle to be be- 
low 10~^eV , with the exception of energy barriers 
determination where the thresholds were 10~^ eV for 
electrons and 10~^ eV for ions. We have checked that 
both setups give the same results within a meV accu- 
racy. 

The slab-supercell considered has been carefully 
tested and a 20 A vacuum along the c axis has been 
adopted to ensure no reciprocal interaction between 
periodical images. We find that using the above set- 
tings the interaction between two adjacent graphite 
layers is ~2 meV, largely within the intrinsic DFT er- 
ror. This result is in agreement with literature data 



Palacios et al. 2008 Yazyev 20081. Indeed, they all 



(Hasegawa and Nishidate 2004 Rydberg et al. 20031 



share the disappearance of one or more carbon p or- 
bitals from the tt — tt* band system, a fact which may 
lead to the appearance of magnetic textures and in- 
troduce site-specific dependence on the chemical prop- 



For this reason a single graphene sheet can also model 
the Bernal (0001) graphite surface, at least as long as 
chemical interactions are of concern. 

The cell size on the surface plane is a fundamental 
parameter for these calculations, since we have found 



3 



U.11IL L.cll 


ft / MT 


H , /A 


P t / pV 
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this work others 


2x2 


0.125 


0.36 0.36^ 


0.75 0.67^ 


3x3 


0.062 


0.42 0.41^ 


0.77 0.76^ 


cluster 


0.045 


0.57^ 


0.76=* 


4x4 


0.031 


0.48 


0.79 0.76^*, 0.85^5 


5x5 


0.020 


0.59 


0.84 0.71^ 0.82'^ 


8x8 


0.008 




0.87* 



Table I Chemisorption energy (Echem) and equilibrium 
height of the C atom above the surface (dp„cfe) for H ad- 
sorption on top of a C atom, for a number of surface unit 
cells, corresp onding to differ ent coverages d. Ref. 1 |Sha| 
and Jackson] p002 ), Ref. 2 [Kerwin and Jackson] (|2008||, 



Ref. 3 Ferro et al. (20 03|, R ef. 4 Duploc k et al 



Ref. 5 |Hornekaer et al.|(j20b6b| ), R ef. 6, Roman et al.|(l2007 l 



Ref. 7 Chen et al. (20071, Ref. 8 Lethinen et al. (2004 I 



(20041 



that chemisorption energies strongly depend on the 
coverage (see below). We choose to use a reasonably 
large 5x5 cell in order to get some tens of meV accu- 
racy while keeping the computational cost as low as 
possible. Even with this size, however, the possibility 
of interactions between images has always to be taken 
into account when rationalizing the data. 



III. RESULTS 

A. Single atom adsorption 

Chemisorption of single H atoms on graphite has 
long been studied since the wor ks of [Jeloaica and Sidis] 
(1999f and Sha and Jackson ( 2002 1 , who first pre- 
dicted surface reconstruction upon sticking. Such a 
reconstruction, i.e. the puckering of the carbon atom 
beneath the adsorbed hydrogen atom, occurs as a con- 
sequence of sp^ — sp^ rehybridization of the valence C 
orbitals needed to form the a CH bond. Since this 
electronic/nuclear rearrangement causes the appear- 
ance of an energy barrier ^--^0.2 eV high, sticking of 
hydrogen atoms turns out to be a thermally activated 
process w hich hardly occurs at and b elow room tem- 
perature (Kerwin and Jackson 20081. 



As already said in Section II, we have re-considered 
adsorption of single hydrogen atoms for different sizes 
of the surface unit cell. We have found that both the 
binding energy and the puckering height are strongly 
affected by the size of the unit cell (see Table |l]| , and 
even the results of the 5x5 cell turns out to be in er- 
ror of about ~30 meV with respect to the isolated 
atom limit estimated by the calculation at 0.008 ML 
coverage (Lethinen et al. 20041. In particular, we 



have found that some cautions is needed in compar- 
ing the height of the carbon atom involved in the bond 
since constraining the neighboring carbon atoms in ge- 
ometry optimization may lead to considerable surface 
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Figure 1 Top panel: total density of states for graphene. 
Bottom panel: density of states for spin-up (positive val- 
ues) and spin-down (negative values) components in a 5x5 
H layer on graphene. 



strain. 

Despite this, we have consistently used the 5x5 cell 
in studying multiple adsorption of hydrogen atoms. 
Indeed, this size allowed us to investigate a number 
of stable configurations involving two, three and four 
adsorbed H atoms, along with the barrier to their for- 
mation, with the same set-up described in Section II. 
Interactions between images do indeed occur for some 
configurations but, as we show below, this does not 
prevent us to get a clear picture of the adsorption 
processes we are interested in. 

In agreement with previous studies we find that hy- 
drogen adsorption can only occur if the substrate is 
allowed to relax. Without relaxation the adsorption 
curves on different surface sites are repulsive, and only 
a metastable minimum is found for the atop position 
( Sha and Jackson 2002 I . Surface relaxation requires 



about 0.8 — 0.9 eV and results in the outward motion 
of the carbon atom forming the CH bond (see Ta- 
blejljand Jeloaica and Sidis (19991; Sha and Jackson 
(2002|). 



In addition, we have investigated the electronic 
substrate properties of the resulting hydrogenated 
graphene, in order to get hints for understanding the 
adsorption process of additional atoms. In Fig. Jl] 
we show the Density of States (DOS) of the 5x5 H- 
graphene equilibrium structure (bottom panel), com- 
pared to that of clean graphene (top panel). It is evi- 
dent from the figure that hydrogen adsorption causes 
the appearance of a double peak in the DOS, sym- 
metrically placed around the Fermi level. This is in 
agreement with rigorous results that can be obtained 
in tight binding theory for bipartite lattices. Indeed, 
Inui et al. ( 1994 1 have shown that for a bipartite lat- 



4 



tice with A lattice sites and ub B lattice sites a suf- 
ficient condition for the existence of mid-gap states is a 
lattice imbalance {ua ^ ub)- In particular, there exist 
ni = \nA — nB\ mid- gap states with vanishing wave- 
function on the minority lattice sites. In H-graphene a 
lattice imbalance results as a consequence of the bond 
with the H atom which makes one of p orbitals no 
longer available for taking part to the tt — tt* band sys- 
tem. There is one mid-gap state for each spin species, 
and the degeneracy is lifted if exchange-correlation ef- 
fects are taken into account, as shown in Figjljfor our 
DFT results. This state has been mapped out in Figj2] 
(left panel) , where we report a contour map of the spin 
density at a constant height 0.47 A above the surface. 
It is clear from the figure that if adsorption occurs 
on a A lattice site the spin-density (due to mainly to 
the above mid- gap state) localizes on B lattice sites. 
The latter now contain most of the 1 ^b magnetiza- 
tion {hb =Bohr magneton) previously carried by the 
H atom species, and a slight spin-down excess on A 
sites results as a consequence of the spin-polarization 
of the lower lying states. This is made clearer in the 
right panel of Fig. [2] where we report the spin-density 
at the same height above the surface as before along a 
rectilinear path joining a number of C atom sites away 
from the adsorption site (see Figjs] for the labels). 
Note that the spin-density decays only slowly with 
the distance from the adsorption site, in agreement 
with theoretical results that suggest that in the case 
of two dimensional graphene this decay corresponds to 
a non-normalizable state with a l/r tail (in contrast 
to non-zero gap substrates such as armchair nanorib- 



bons where mid-gap states are normalizable) ( Pereira 



et al. 20061. With our unit cell the effect of the in- 
teraction with the images is already evident at rather 
short distances, but as we show below, this effect has 
no infiuence on the interpretation of the results. Note 
also that this spin pattern is common to other 'de- 
fects' (e.g. vacancies, voids and edges) which have 
been known for some time to strongly modify the elec- 
tronic properties of graphene and graphene-like struc- 



magnetic structures (Pereira et al. |2006, ,Yazyev and 


Helm| |2007| |Pisani et al.l 2008| |Pereira et al.| 2008| 


Palacios et al.| |2008| Yazyev 


'20C 
, Rr 


i8l Lethinen et al.' 


2004| Mizes and Foster 1989 


iffieux et al. 




2000. 


Kusakabe and Maruyama| 2003 


Jiang et al. 


2007 


Yazyev et al. 


2008 


Yazyev and Katsnelson 


20081. In 



particular, in a recent, comprehensive study, Palacios 
et al., using a mean-field Hubbard model for graphene, 
have clarified the appearance of magnetic textures as- 
sociated to vacancies and predicted the emergence of 
magnetic order (Palacios et al. 20081. Their model 



also suits well to 'defects' such as the presence of the 
adsorbed hydrogen atoms. 

From a chemical point of view the above spin pat- 
tern (and the resulting magnetic properties) arise from 
the 'spin-alternation' typical of tt conjugated com- 





Figure 2 Spin density 0.47 Aabove the graphene surface 
after adsorption of a hydrogen atom. Left: contour map 
with red/blue lines for spin-up/spin-down excess respec- 
tively. Right: spin-density at the same height as on the 
left panel, along a path joining the C atoms (for the labels 
see Figj5|. 



pounds. This behavior is easily understood in terms 
of resonant chemical structures, such as those shown 
m Figj3] for a coronene molecule. In this and analo- 
gous Polycyclic Aromatic Hydrocarbons (PAHs), the 
TT electron system can be described as a resonant 
combination of conventional, alternated double bond 
structures, like the one shown in the upper panel of 
Figjs] (see a). Once a hydrogen atom has been ad- 
sorbed on the surface, an unpaired electron is left 
on one of the neighboring C atoms (6, left panel), 
which can subsequently move in each of the carbon 
atoms belonging to a sublattice V3xv/3R30° by 'bond- 
switching' (see h,c). Spin-alternation arises from the 
'resonant' behavior of an unpaired electron in a po- 
sition (the nearest neighbor one) with respect to a 
double bond: such 'resonance' can be naively viewed 
as the spin re-coupling of the unpaired electron with 
the electron on the neighboring site, a process which 
sets free a second electron on the same sublattice. 

This picture, despite its embarrassing simplicity, 
can be put on firm grounds in the context of the Va- 
lence Bond (VB) theory of chemical bonding (see e.g. 



Raimondi et al. (19851; Cooper et al. 


(19871; 


Gerratt 


et al. 


( 1997} ; Li and McWeeny (|2002 1 ; 


Cooper 


(20021). 



Focusing on the tt electron system, this can be done 
with the help of a simple (correlated) VB ansatz for 
the N electron wavefunction of the tt cloud, namely 



'^SN = -4(0102. .(/)Ar6sJv) 



(1) 



where A is the antisymmetric projector, 0^ = 0i('") 
for i = 1, TV are (spatial) orbitals accommodating the 
TV electrons, and Qsn is a N electron spin-function 
with spin quantum number S. The latter is usually 
variationally optimized by expansion on a spin func- 
tion basis, 



e 



SN 



E 



SN:k 
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Figure 3 (a) The vr resonating chemical model for a 
grapheme surrogate (coronene). (b), (c) Spin-alternation 
after hydrogen adsorption. 



where is the dimension of the spin-subspace of 
eigenfunctions of with eigenvalue S{S + 1) and 
given magnetization^. Among these basis-functions 
the 'perfect pairing' set devised by Rumer, though 
non-orthogonal, is chemically appealing since for a 
given S and Mg = S the total magnetization is given 
by 25' electrons coupled at high spin, the remaining 
N — 2S being accommodated in {N — 2S) /2 singlet- 
coupled pairs (see Simonetta et al. ( 1968| ). Then, if 
the orbitals are localized on the atoms, the resulting 
wavefunction 



CkA{<pi(l)2-(f>NQsN;k) = ^ Ck^ 



SN;k 



is a superposition of conventional 'structures' 'ifsN-k 
describing pairs of atom-centered, singlet-coupled or- 
bitals (i.e. Lewis chemical bonds and lone pairs) and 
unpaired electrons. 'Classical' molecules require just 
one perfect-pairing spin function coupling those pairs 
of orbitals with substantial overlap. Less conventional 
molecules, such as tt conjugated systems, need a true 
superposition of two or more spin structures, since the 
energy gain (also known as resonance energy) in al- 
lowing such superposition is particularly important in 
these cases. Correspondingly, the classical Lewis pic- 
ture of chemical bonds is extended to account for the 
resonance phenomenon, as shown in Fig. |3]with dou- 
ble ended arrows indicating superposition of chemical 
structures. 



1 is given by the expression = Ar!/(7V/2-|-5-l-l)!/(iV/2- 
Sy.(2S + 1) and does not depend on the value Ms of the 
spin-projection zS along the axis z, since these subspaces are 
isomorphic to each other. 



Early applications of the theory, starting from the 
landmark work of Heitler and London, used frozen 
atomic orbitals. In modern, ab initio use of the the- 
ory both the spin-coupling coefficients Ck and the or- 
bitals can be variationally optimized, even when us- 
ing a number of configurations in place of the single 
orbital product appearing in eq.Q (see e.g. Li and 
McWeeny (2002 l), in close analogy to what is done in 
molecular orbital theory with the MultiConfiguration 
Self-Consistent Field (MCSCF) approach. The inter- 
esting thing is that these optimized orbitals, as a con- 
sequence of electron correlation, are usually (if not al- 
ways) localized on atomic centers and are only slightly 
polarized by the environment ( Cooper et al.[ 1987 



Gerratt et~aL] |1997[ [Cooper et al.[ 1 1 99 1| , thereby 



supporting the interpretation of the simple wavefunc- 
tion of eq. ([ij as a quantum-mechanical translation of 
Lewis theory of chemical bond. This is true, in par- 
ticular, for the benzene molecule, the prototypical tt 
resonant system, where six, p-like orbitals are mostly 



coupled by two, so-called Kekule structures (Cooper 
et al l [T9§7l [Tantardini et al.| [19771 [Cooper et "aL 
1986P. 



From a physical point of view, wavefunction ([T| 
generalizes to N electron systems the Heitler-London 
ansatz forming the basis for the Heisenberg model 
of magnetism in insulators. In addition, if the or- 
bitals are allowed to be 'polarized', band-like behav- 
ior can be accommodated, along with collective spin 



excitation, as in the Hubbard model (Hubbard 19631 



which has been finding widespread use in investigating 
graphitic compounds. The fact that Hubbard model, 
and its Heisenberg limit, can be derived by suitable 
approximations to Valence-Bond ansatz has long been 
known in the chemical literature, esp ecially in connec 
tion to TT resonant systems (see e.g. Wu et al. (2002 
2003 1 and references therein. Hubbard model is also 



known as Pariser-Parr-Pople model in the chemical 
literature 



after Pariser and Parr (Pariser and Parr 
1953a|b[ | and after Pople ( [Pople| ]T953D 7 



We can 

roughly say that Heisenberg models correspond to the 
'classical' valence theory developed by Heitler, Lon- 
don, Pauling and Van Vleck in the twenties which put 
the basis for explaining chemical bond using frozen 
atomic orbitals, whereas Hubbard models arise from 
the modern version of theory, started with Coulson 
and Fisher and pushed forward by Gerratt and oth- 
ers, who used 'polarizable' orbitals in the same spin 
scheme set up in the original theory (Cooper et al. 



2 For 5 = and A'" = 6 the set of five linearly independent 
Rumer structures is given by two Kekule structures and the 
three additional 'Dewar' structures. A resonance energy ~0.8 
eV can be computed when using two Kekule structures in 
place of one, whereas only some tenths of meV are gained 
when full optimiz ation of the spin function is performed ( |Bon-| 
fanti et al. 1 120081. 
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1987 19861. 



In the following Sections, we will use the above 
wavefunction (eq.Q) as a simple guide to interpret 
the results of our first -principles calculations, keeping 
in mind its connections with the traditional chemical 
picture on the one hand and the Hubbard model on 
the other. As we will show in the following, even at 
this qualitative level, a number of useful insights can 
be gained from such a picture. As a first example we 
can reconsider adsorption of a single hydrogen atom 
on graphene. In a diabatic picture (i.e. when con- 
straining the spin-coupling to the Kekule structures of 
Figjs] panel a) the interaction between graphene and 
the incoming H atom is expected to be repulsive, since 
no electron is available to form the CH bond. On the 
other hand, a low lying spin-excited state correspond- 
ing to a Dewar-like structure (which has two, singlet- 
paired electrons on opposite, no-overlapping end of a 
benzene ring) would give rise to an attractive, barri- 
erless interaction. At short-range, then, an avoided 
crossing between the two doublet curves occurs which 
signals the spin-transition leading to bond formation, 
even though this can lead only to a metastable state if 
surface reconstruction is not allowed, as indeed found 



in DFT calculations (e.g. see Fig. 2 in (Sha and Jack 



2002 1). Actually, in this case the situation is a 



bit more complicated since a slightly lower-lying state 
in the triplet manifold (obtained by spin-flipping the 
above spin-excited Dewar-like structure) contributes 
to the same doublet manifold. Valence Bond calcula- 
tions on the simpler benzene- H system confirms this 
picture ([Bonfanti et al.||2008|, see FigjI] 




Figure 4 Interpretation of the sticking barrier as an 
avoided crossing between chemical structures. Valence 
bond results for the benzene-H system, from Ref. (Bon- 



fanti et al. 20081. Solid black and red circles for the 
ground [CeHeC Aig) + H(^S)) and the first excited states 
{CeHei^ Biu) + H(^S)), as obtained at the single-orbital- 
string level of eqjl] with orbital optimization. Quasi- 
diabatic results are obtained by properly constraining the 
spin space: Kelule structures only (lower right and up- 
per left insets) for empty black circles; structures in the 
lower left inset for empty red circles. Also shown in the 
upper right inset the main (Dewar-like) structures needed 
to described the ^Biu state of benzene. 



B. Secondary adsorption 

Next we consider adsorption of a second atom on 
the different sites A(n) {n — 1,6) and B(n) {n — 1,6) 
shown in FigjS] with a first adsorbed H atom on site 
A(0). For each site we have investigated the ground 
spin manifold by allowing full relaxation of the magne- 
tization. In addition, in most of the cases, we have also 
performed magnetization-constrained calculations in 
order to get insights on both the singlet and the triplet 
states arising from the interaction between the dou- 
blet H-graphene ground-state and the second H atom 

The results for the binding energies are reported in 
Table [nl along with the site-integrated magnetizations 
(Ms/) and the total magnetization M. Site-integrated 
spin-densities have been obtained by integrating the 
spin-density on a small cylinder (of radius half of the 
C-C distance in the lattice) centered on each site, and 
can be considered a rough measure of the total spin ex- 
cess available on the site. This quantity behaves very 
similar to the spin-density itself, decreasing in magni- 
tude when increasing the distance from the adsorption 
site, separately for each sublattice. Some exceptions 



Position 


Msi 1 MB 


^bind / eV 


M/ms 


^binH / eV 


B(l) 


0.109 


1.934 





0.933 


A(l) 


-0.019 


0.802 


2 


0.575 


B(2) 


0.085 


1.894 





0.828 


A (2) 


-0.017 


0.749 


2 


0.531 


B(3) 


0.040 


1.338 





0.646 


A (3) 


-0.016 


0.747 


2 


0.570 


B(4) 


0.076 


1.674 







A (4) 


-0.016 


0.747 


2 


0.573 


B(5) 


0.023 


1.033 





0.590 


A (5) 


-0.014 


0.749 


2 


0.531 


B(6) 


0.028 


1.110 





0.545 


A (6) 


-0.015 


0.787 


2 





Table II Binding energies (^^^ind) '^'^^ secondary adsorp- 
tion to form the H-pairs shown in FigjS] along with the 
site-integrated magnetizations (Mg/) before adsorption, 
and the total ground-state magnetization (M) after ad- 
sorption obtained when fully relaxing the magnetization. 
Also reported the binding energies obtained when the mag- 
netization is constrained to M = 0, 2 /is for A and B sites, 
respectively. See text for details. 
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Figure 5 The graphene unit cell used for the calculations 
with A (blue) and B (red) lattice sites indicated. Also 
indicated is the path used for Figj2] ^(0) being the first H 
adsorption site. 



are worth noticing, namely the A(0)-B(5) pair, and 
are due to the cumulative effect of next-neighbors im- 
ages. Notice, however, that despite their possible arti- 
ficial nature, results corresponding to any lattice sites 
when viewed as a function of the site-integrated mag- 
netization give insights into the adsorption process. 

A quick look at the Table |ll] reveals that the two 
sublattices A and B behave very differently from each 
other, as the spin-coupling picture of Fig. |3] (panels 
6,c) suggests. Roughly speaking, adsorption on B lat- 
tice is preferred over that on the A lattice. The bind- 
ing energies are much larger than the first adsorption 
energy reported in Table |l] (they can be as large as 
twice the adsorption of the first atom) , and give rise 
to a final unmagnetized state. In contrast, the binding 
energy for adsorption on a A lattice site is comparable 
to that of single-H adsorption, and the ground-state 
of the H-pair on graphene is a triplet (M=2 /is). 



These findings agree with Lieb's theorem (Lieb 



1989 1 for the repulsive Hubbard model of a bipartite 
lattice and a half-filled band, which states that the 
ground-state of the system has S =l/2\nA — tib]- In 
such model, the electronic state of the system would 
be described hy N — 2 p orbitals {N being the original 
number of sites), and riB — riA — N/2 — 1 if adsorp- 



> 

0) 



LU 



3.00 



2.00 



1.00 



0.00 




0.00 



0.05 0.10 



0.15 



B 



Figure 6 Binding energies for secondary H adsorption as 
a function of the site-integrated magnetization, for sin- 
glet (red squares) and triplet (blue squares) states. Black 
square is the data point for single H adsorption. Also 
shown a linear fit to the data set (solid line) and the H 
binding energy to form some 4-atom clusters from 3-atom 
ones (red and blue circles for final singlet and triplet states, 
respectively). See text for details. 



tion of the second hydrogen atom proceeds on the B 
lattice (to form what we can call AB dimers), whereas 
riB = nA + 2 — N/2 if it proceeds on the A lattice (to 
form A2 dimers). The results are also consistent with 
the VB framework sketched in Subsection IIII.AI with 
reference to Fig. |3] (panels b,c), it is clear that when 
a H atom adsorbs on an B site its electron readily 
couples with the unpaired electron available on the B 
sublattice, whereas when adsorption occurs on an A 
site two electrons are left in excess on the B sublattice, 
and they more favorably couple at high spin. 

The relationship between the available unpaired 
electron density at a given site and the binding energy 
of adsorbing a second H atom can be made clearer by 
reporting the energy data of Table |ll] as a function 
of M5/. This is shown in Figj6] for both the singlet 
and triplet states of the dimers, along with the value 
for the first H adsorption (data point at M5'/=0). 
It is clear from the figure that, with the exception 
of the value for the orifto-dimer (A(O)S(l) in Figj5] 
This value has been excluded from the linear regres- 
sion shown in Fig|6|, a linear relationship between the 
binding energy and the site integrated magnetization 
well describes the situation, and the binding energy 
for single H adsorption fits well to this picture. 

This is again consistent with the chemical model, as 
long as the site-integrated magnetization is a measure 
of the unpaired electron density available. Accord- 
ing to Section |III.A| adsorption of the first hydrogen 
atom arises from the energy balance between a 'local- 
ization energy' (the spin excitation needed to set free 
an unpaired electron on the given lattice site), the 
spin-pairing forming the bond, and the surface recon- 
struction energy. The same is true for adsorption of a 
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0.15 



Figure 7 Barrier energies for secondary atom adsorption 
on the B lattice sites as a function of the site-integrated 
magnetization. Linear regression of the data omits the 
values for forming para and ortho pairs (two rightmost 
points in the graph). 



second atom: localization energy takes only a slightly 
different form than before because an unpaired elec- 
tron is already available in one of the two sublattices^ , 
but surface reconstruction energy is not expected to 
depend on the adsorption site. Then, adsorption en- 
ergies depend on the electronic properties only, and 
the linear behavior observed for singlet-state dimers 
in Fig. |6] suggests that the energy needed to localize 
the unpaired electron on a given site decreases linearly 
when increasing the unpaired electron density avail- 
able. Notice that negative values of M5/ (as found at 
A sites), correspond to a spin excess parallel to that of 
the incoming H electron, and for these sites localiza- 
tion of an unpaired electron with an antiparallel spin 
requires increasingly more energy when the (magni- 
tude) of the spin-density increases, since this can only 
be achieved by adding one electron to the site. On 
the other hand, when a triplet dimer is formed upon 
adsorption the H electron does not make use of the 
unpaired electron available, and adsorption energies 
are all around ~0.8 eV, i.e. of the order of the first 
H adsorption. The effect of surface relaxation is only 
seen in forming the ortho-dimei, where few tenths of 
eV more than the single H relaxation energy are re- 
quired because of the closeness of the two hydrogen 
atoms. 

Analogous linear behavior can be found when con- 



In terms of the wavefunction of eq.([T| this localization energy 
can be defined by observing that the structures in which the 
spin-up density localizes on the {N — 1) — th site (N even) 
correspond to 't = A((/)i02--</'JV-i0;^J^)i where ©^"'^ is 
constrained to have the form ©j'^J^ = {c^Q^ +C2@'^ Q^)a, 
whereas the ground-state spin function comprises additional 
contributions from ©^j"^/3 structures. 



sidering the computed energy barrier to sticking as 
a function of the site-integrated magnetization, as 
shown in Fig. [7] for AB dimers. This agrees with the 
above localization energy and with the common ten- 
dency for a linear relationship between the binding 
and the barrier energies for activated chemical reac- 
tions (Br0nsted-Evans-Polayni rule). Exceptions are 
given by the ortho- dimer considered above and by 
the para- dimer {A{Q)B{2)). The latter, in particular, 
shows no barrier to adsorption, in agreement to pre- 
vious theoretical works, and this fact forms the ba- 
sis of the so-called preferential sticking mechanism. 
This mechanism was first suggested by Hoernaker et 
al. (iHornefcer et all |2006bf who looked at the STM 



images formed by exposing Highly Oriented Pyrrolitic 
Graphite (HOPG) samples to a high-energy (1600- 
2000 K) H atom beam and observed formation of 
stable pairs, confirmed by first -principles calculations 
(Hornekaer et al. 2006b I. Our results suggest that 



barrierless adsorption on the para site is a consequence 
of both favorable electronic and nuclear factors. 

We therefore find that formation of AB dimers is 
both thermodynamically and kinetically favoured over 
formation of A2 dimers and single atom adsorption. 
This agrees with current experimental observations 
which show evidence for clustering of hydrogen atoms 
at all but very low (< 1%) coverage conditions. In 
addition, we notice that the dimers identified so far 
Hornekffir et al. 2006a|b I are all 



(Andree et al. 2006 



of the AB type. 



C. Further adsorptions 

We consider in this section results concerning for- 
mation of cluster of three and four atoms. In these 
cases, the number of possible configurations is quite 
large and therefore we limit our analysis to a few im- 
portant cases. Following analogous notation recently 
introduced for defects by Palacios et al. (20081, we use 
the 'chemical formula' AnB^ to denote a cluster with 
n H atoms in the A lattice and m H atoms in the B 
lattice. According to Lieb's theorem and to the tt res- 
onance picture, we expect that the ground electronic 
state has |n — m| unpaired electrons. We have consid- 
ered a number of A2B2, A2B, A3B1 and A3 clusters, 
and found that their ground-state has 0, 1, 2 and 3 
liB of magnetization, respectively, in agreement with 
the expectation. 

Three atom clusters have been obtained by adding 
one hydrogen atom either to a para dimer or to a meta 
dimer, i.e. A{0)B{2) and A{0)A{1) with the labels of 
Figj5] respectively. The binding energies of a third hy- 
drogen atom to a para dimer structure are reported in 
Tab jIII[ since they all are of A2B type, the total mag- 
netization for the resulting structures is 1 /i^. A look 
at Tab jTn] reveals that adsorption to a third hydrogen 
atom parallels that of the first H. This is consistent 
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Position 




A(2) 


1.516 


B(3) 


0.847 


A(3) 


0.727 


B(4) (=A(5)) 


0.971 


A(4) (=B(6)) 


0.821 


B(7) 


0.727 


B(8) 


1.301 



Table III Binding energies (^JbiiKj) '^'^^ addition of a third 
H atom to the para dimer structure y4(0)i3(2) on the sites 
indicated in the first column (labels from Fig. [5|. 




Figure 8 Spin-density 0.40 A above the surface for two 
three-atom clusters. Contour map with red/blue lines 
for spin-up/spin-down excess respectively. Left and right 
panel for an A2B and a A3 cluster, respectively. 



with the TT resonance picture, since AB dimers do not 
have unpaired electrons, and therefore show no pref- 
erence towards any specific sublattice position. There 
are of course exceptions, notably the values for ad- 
sorption onto A(2) and i?(8) lattice sites, and these 
can be reasonably ascribed to the effect of surface 
relaxation. Indeed, relaxation energies per atom in 
'compact' clusters may considerably differ from the 
value of the single H atom, being always of the order 
of the binding energies themselves (~ 0.8 eV). Simi- 
lar conclusions hold when adding a third H atom to 
the (magnetic) meta dimer A(0)A(1): adsorption on 
B lattice sites is strongly favored {Ehind = 1-2 — 1.9 
eV) and produces doublet structures (Af = \\xb), 
whereas H atoms bind to A lattice sites with an energy 
~ 0.7 — 0.8 eV and produce highly magnetic structures 
(M = 3 [Ib)- Energy barriers to adsorption follow the 
same trend: preliminary calculations show that, with 
few exceptions, barriers to sticking a third H atom 
compare rather well with that for single H atom ad- 
sorption for the processes AB A2B and A2 A3, 
and may be considerably smaller for A2 ^ A2B ones. 

In addition, again consistently with tt resonance pic- 
ture, we found that all the considered 3-atom struc- 
tures, with one or two unpaired electrons, show an 
alternation pattern in their spin-density maps. As 
an example, FigjS] reports the spin-density maps for 





Mqt I U 
01 / t^rS 


Eu- J / eV 
'^bind / 


M //; H 


B(9) 


-0.0180 


1.103 


2 


A{7) 


0.0471 


1.331 





B(6) 


-0.0151 


0.727 


2 


A(8) 


0.0325 


1.210 





5(10) 


-0.0134 


0.723 


2 


A{9) 


0.0326 


1.201 






Table IV Binding (E^ind) energies for adsorption to form 
H-quadruples from the A{0)B{2)B{8) cluster, along with 
the site-integrated magnetizations (Ms/) and the total 
ground-state magnetization (M), before and after adsorp- 
tion, respectively. See FigjS] for atom labels. 

an A2B (left panel) and an A3 (right panel) cluster. 
Analogously to Subsection |III.B| we find that anal- 
ysis of these spin-density maps gives insights to the 
adsorption properties of a fourth hydrogen atom. Ta- 
ble |IV| for example, reports binding energies to form 
some 4-atom clusters from the stable A{0)B{2)B{8) 
one, the final total magnetization of the resulting 
structures and the values of the corresponding site- 
integrated magnetization before adsorption. The com- 
puted binding energies compare rather well with the 
dimer values, as can be seen in Fig. [6] where it is 
clear that the results fit well to the same linear trend 
obtained before. Few exceptions are for compact clus- 
ters where substrate relaxation does play some role. 
With such exceptions in mind, our results suggest that 
adsorption of hydrogen atoms on magnetic graphitic 
substrates (such as those obtained by adsorbing an 
odd number of H atoms), for a given final spin-state, 
depends on the local spin-density only. 



IV. SUMMARY AND CONCLUSIONS 

In this work we have presented results of extensive 
first -principles calculations of the adsorption proper- 
ties of hydrogen atoms on graphite. A number of pos- 
sible configurations involving one, two, three and four 
atoms on the surface have been considered and bar- 
rier energies have been computed for some of them. 
We have found that adsorption of hydrogen atoms 
is strongly related with substrate electronic proper- 
ties, and used the chemical model of planar tt conju- 
gated systems to rationalize the data. The connec- 
tion between this model and the valence theory of 
chemical bond on the one hand, and Hubbard mod- 
els on the other hand, has been emphasized in Sec- 
tion |III.A| and used at a qualitative level to rational- 
ize our findings. In this way, one prominent feature 
of defective graphitic substrates, i.e. the possibility 
of forming ordered (microscopic) magnetic patterns, 
turns out to be related to the spin-alternation typi- 
cal of TT resonant systems. We have also invariably 
found in the cases considered that Lieb's theorem for 
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repulsive Hubbard models can be used to predict spin 
alignment in ground-state graphitic structures. 

Adsorption of single H atoms has been known for 
some time to be an activated process, with an energy 
barrier to sticking (~ 0.2 eV) high enough to prevent 
adsorption at ambient conditions. Adsorption of a sec- 
ond atom more favorably occur on the y/Sxy/smO" 
sublattice where spin-density localizes, and may pro- 
ceed without barrier if it occurs on the so-called para 
site. This preferential sticking has been recently sug- 
gested by experimental and theoretical observations 
(Hornekffir et al. ( 2006a|b |). We extended the lat- 
ter analysis by considering a large number of possible 
dimers and found that (i) binding (barrier) energies 
generally increase (decrease) linearly as a function of 
the site-integrated magnetization, and (ii) adsorption 
properties of the ortho and para sites are slightly at 
variance with linear trends, thereby suggesting that 
substrate relaxation plays some role in these cases. 

When considering addition of a third atom we found 
that the adsorption energetics of the incoming H atom 
is similar to that of the first one (i.e. a barrier ~ 0.2 
eV high and a chemisorption well ~ 0.8 eV deep), un- 
less we start with a 'magnetic' dimer in which the two 
atoms are adsorbed in the same sublattice. (These 
structures, however, are kinetically and thermody- 
namically unfavored with respect to the unmagnetized 
AB configurations). This is in agreement with the 
chemical model, which predicts an open-shell configu- 
ration for A2 dimers and a closed-shell one (with par- 
tial restoring of the tt aromaticity) for AB ones. These 
results, therefore, suggest that preferential sticking 
alone cannot provide any catalytic route to molecu- 
lar hydrogen formation on graphite. 

Finally, we have considered adsorption energetics in 
forming clusters of four atoms, and re-gained the same 
picture obtained in forming pairs, namely that adsorp- 
tion is strongly biased towards the sublattice in which 
the spin-density localizes. Actually, the resulting en- 
ergetics fits well to the linear behavior with respect 
to the site-integrated magnetization already found for 
dimer formation. Such a linear relationship suggests 
that the energy needed to localize the unpaired elec- 
tron on a given lattice site decreases linearly when 
increasing the site-integrated magnetization, at least 
in the range of values covered by this study. Inter- 
estingly, this behavior suggests that if we were able 
to tune the magnetization of the substrate we could 
control the adsorption dynamics of H atoms. 

Overall our results, consistently with the tt res- 
onance picture, suggest that the thermodynamically 
and kinetically favored structures are those that mini- 
mize sublattice imbalance, i.e. those AnB„i structures 
for which ni — \n—m\ is minimum. The latter number 
71/ is also the number of mid-gap states in single par- 
ticle spectra which, accordin g to the Hu nd-like rule 
provided by Lieb's theorem (Lieb 19891, is directly 
related to the total spin of the system, S = n//2. 



which is therefore at minimum in the favored struc- 
tures. Notice that however small the S value can be, 
this result does not preclude the existence of local 
magnetic structures, antiferromagnetically coupled to 
each other. The case of an AB dimer with two atoms 
very far from each other provides such an example. 
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